
install.packages("mediation")
library("mediation")


### Table 6 - Model 1 ###

newdata <- RawData[ which(RawData$ConditionCode !=3),]
newdata2 <- newdata[ which(newdata$ConditionCode !=4),]
newdata3 <- newdata2[ which(newdata2$ManiCheck1 == 1 | newdata2$ManiCheck2 == 2 | newdata2$ManiCheck1 == 0),]

model.m <- lm(STAXI_Mean ~ ConditionCodeY + DummyCountryCodeUS + DummyCountryCodeUK + PoliticalPosition + Year_4 + Gender, data = newdata3)
model.y <- lm(Retaliation_Reduced_Mean ~ STAXI_Mean + ConditionCodeY + DummyCountryCodeUS + DummyCountryCodeUK + PoliticalPosition + Year_4 + Gender, data = newdata3)
out.1 <- mediate(model.m, model.y, sims = 1000, boot = TRUE, treat = "ConditionCodeY", mediator = "STAXI_Mean")

summary(out.1)
plot(out.1)



### Sensitivity Analysis - Supplementary Analyses - Online Appendix F - Left Model ###

med.cont <- mediate(model.m, model.y, sims=1000, treat = "ConditionCodeY", mediator = "STAXI_Mean")
sens.cont <- medsens(med.cont, rho.by = 0.01)

summary(sens.cont)

plot(sens.cont, sens.par = "rho")
plot(sens.cont, sens.par = "rho", xlim = c(-.5,.5), ylim = c(-.5,.5))








### Table 6 - Model 2 ###

Newdata9 <- RawData[ which(RawData$ConditionCode !=1),]
newdata8 <- Newdata9[ which(Newdata9$ConditionCode !=2),]
newdata7 <- newdata8[ which(newdata8$ManiCheck1 == 1 | newdata8$ManiCheck2 == 2 | newdata8$ManiCheck1 == 0),]


model.m2 <- lm(STAXI_Mean ~ ConditionCodeZ + DummyCountryCodeUS + DummyCountryCodeUK + PoliticalPosition + Year_4 + Gender, data = newdata7)
model.y2 <- lm(Retaliation_Reduced_Mean ~ STAXI_Mean + ConditionCodeZ + DummyCountryCodeUS + DummyCountryCodeUK + PoliticalPosition + Year_4 + Gender, data = newdata7)
out.2 <- mediate(model.m2, model.y2, sims = 1000, boot = TRUE, treat = "ConditionCodeZ", mediator = "STAXI_Mean")

summary(out.2)
plot(out.2)



### Sensitivity Analysis - Supplementary Analyses - Online Appendix F - Right Model ###

med.cont2 <- mediate(model.m2, model.y2, sims=1000, treat = "ConditionCodeZ", mediator = "STAXI_Mean")
sens.cont2 <- medsens(med.cont2, rho.by = 0.01)

summary(sens.cont2)

plot(sens.cont2, sens.par = "rho")
plot(sens.cont2, sens.par = "rho", xlim = c(-.5,.5), ylim = c(-.5,.5))



